Strongly nonlinear waves in capillary electrophoresis 
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In capillary electrophoresis, sample ions migrate along a micro-capillary filled with a background 
electrolyte under the influence of an applied electric field. If the sample concentration is sufficiently 
high, the electrical conductivity in the sample zone could differ significantly from the background. 
Under such conditions, the local migration velocity of sample ions becomes concentration dependent 
resulting in a nonlinear wave that exhibits shock like features. If the nonlinearity is weak, the sample 
concentration profile, under certain simplifying assumptions, can be shown to obey Burgers' equation 
(S. Ghosal and Z. Chen Bull. Math. Biol. 2010 72(8), pg. 2047) which has an exact analytical 
solution for arbitrary initial condition. In this paper, we use a numerical method to study the 
problem in the more general case where the sample concentration is not small in comparison to the 
concentration of background ions. In the case of low concentrations, the numerical results agree 
with the weakly nonlinear theory presented earlier, but at high concentrations, the wave evolves in 
a way that is qualitatively different. 

PACS numbers: 87.64.Aa,87.15.Tt 
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I. INTRODUCTION 

In capillary electrophoresis (CE) , separation of charged 
molecular species is accomplished by exploiting the dif- 
ferential migration of ions in a narrow channel (10-100 
//m) in which a strong electric field (~ 100 V/m) is ap- 
plied in the axial direction 0, Q • The sample ions exist 
in solution in an electrolytic buffer which is referred to as 
the background electrolyte (BGE). Separation is accom- 
panied by the competing process of diffusive spreading in 
the axial direction which causes peak dispersion. Disper- 
sion reduces resolution of the separation and may lower 
the peak concentration to below the detection threshold. 
It is therefore detrimental. Any effect that tends to in- 
crease axial spreading over the minimum imposed purely 
by molecular diffusion in the axial direction is referred 
to as "anomalous dispersion" 0] . The transport problem 
of ions in the capillary is of considerable interest as it 
determines the amount of dispersion of the sample peak. 

In this paper we are concerned with an effect known as 
"electromigration dispersion" (EMD) that causes signif- 
icant anomalous dispersion when the ratio of sample to 
background ion concentration becomes large enough. For 
this reason it is also known as the "sample overloading 
effect". In CE, it is desirable to have the sample con- 
centration at the inlet as high as possible (to ensure that 
even trace components are within detectable limits) and 
buffer conductivity as low as possible (to minimize Joule 
heat), so that the limitation imposed by EMD quickly 
becomes significant [4j. 

The physical mechanism of EMD may be explained 
roughly in the following way: when the concentration of 
sample ions is sufficiently high in comparison to that of 
the background electrolyte, the local electrical conduc- 
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tivity of the solution is altered in the region around the 
sample peak. However, charge conservation requires the 
electric current to be the same at all points along the 
axis of the capillary. If diffusion currents due to concen- 
tration inhomogeneities are ignored for the moment, it 
follows, that the electric field must change axially. This 
is because Ohm's law, taken together with current con- 
servation, implies that the product of the conductivity 
and electric field must remain constant along the cap- 
illary. The axially varying electric field then alters the 
effective migration speed of the sample ions, which in 
turn alters its concentration distribution. Thus, in the 
continuum limit, the concentration of sample ions is de- 
scribed by a nonlinear transport equation. As expected, 
the CE signal exhibits features reminiscent of nonlinear 
waves familiar from other physical contexts p| H[ . 

A one dimensional nonlinear hyperbolic equation for 
the sample ion concentration may be derived using sim- 
plifications that arise from assuming local electroneutral- 
ity and from neglecting the diffusivity of ions [7fl . The re- 
striction to zero ionic diffusivities was recently removed 
by Ghosal and Chen They considered the minimal 
model of a three ion system - the sample ion, a co-ion and 
counter-ion. The diffusivities of the three ionic species 
were assumed equal, though not necessarily zero. The 
sample ion concentration was then shown to obey a one 
dimensional nonlinear advection-diffusion equation which 
reduced to Burgers' equation if the sample concentration 
was not too high relative to that of the background ions. 

In this paper, we focus on the minimal three ion sys- 
tem considered by Ghosal and Chen [f| but we do not 
assume that the concentration of sample ions is small. 
Local electro-neutrality is however an excellent approxi- 
mation in CE systems, since characteristic length scales 
are much larger than the Debye length which is on the 
order of nanometers. We therefore exploit it to reduce 
the numerical stiffness of the coupled ion transport equa- 
tions. We identify a small number of parameters that 
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primarily determine the system evolution and study the 
dynamics for a representative range of these parameters. 
We show that at low concentrations the peak evolves in 
accordance with the weakly nonlinear theory [sj], but at 
high enough concentrations, the dynamics of peak evolu- 
tion is qualitatively different as the system is dominated 
by the nonlinearity. Surprisingly, in the strongly nonlin- 
ear regime, the peak breaks up into two zones marked by 
a critical concentration (<fr = <p c ) and separated by a dif- 
fusive boundary. The high concentration zone {<f> > (f> c ) 
remains quasi-stationary whereas the low concentration 
zone propagates forward forming a "surge front" super- 
ficially resembling nonlinear wave phenomena familiar in 
the context of water waves, such as a river bore Q. The 
critical concentration (</> c ) can be predicted by a simple 
argument based on flux conservation. At late times, dis- 
persion ensures that concentrations throughout the do- 
main get smaller and the peak once again may be de- 
scribed by Burgers' equation. The complex nonlinear 
behavior is a consequence of the nonlinearity inherent 
in the Nerst-Planck equations of ion transport, just as 
the behavior of large amplitude water waves arise from 
the nonlinear nature of the Navier-Stoke's equations of 
hydrodynamics. 



II. NUMERICAL SIMULATIONS 

We set up and solve numerically an idealized problem 
in which a sample peak migrates in a background elec- 
trolyte. The channel walls are assumed charge neutral, 
so that electro-osmotic flow is absent [18j . Further, local 
electro-neutrality is invoked which enables us to express 
the electric field in terms of the instantaneous concentra- 
tion distributions rather than solve the Poisson's equa- 
tion for the electric potential. This considerably simpli- 
fies the numerical work as the Poisson's equation is stiff 
on account of the smallness of the Debye length. Thus, 
the problem is reduced to solving a set of one dimen- 
sional coupled partial differential equations for the ion 
concentration fields. 



(Di/ui = D/u = ksT where ks is Boltzmann's constant 
and T is the absolute temperature). Note however, since 
the valence Zi are different, the electrophoretic mobilities 
of the species fii = Zieu are not identical. Then the only 
parameters in the problem are Pe = vwo/D, which may 
be regarded as a "Peclet number" based on the electromi- 
gration velocity v, and the two valence ratios z n j z, z p /z, 
where z pi z n and z are respectively the valence of cations, 
anions and sample. We present results for two values of 
the Peclet number: Pe = 100 and 200 and fix the va- 
lence ratio at z : z p : z n = 1 : 2 : — 1. For other values 
of these parameters the results are qualitatively similar. 
The parameter of greatest interest is the degree of sam- 
ple loading or the amplitude of the initial peak. The 
shape of the wave is insensitive to initial conditions, so 
for convenience we take the initial peak shape to have a 
rectangular [l9j profile of height c m and width 2wq cen- 
tered at x — 10u>o- This is also the most common initial 
shape encountered in practice where the sample is intro- 
duced by electrokinetic injection. The degree of sample 
loading is conveniently characterized Q in terms of the 
quantity 

/+oo r+oo 
(f)(x,t) dx = — dx, (1) 

- OO «/ — OO ("71 

which has units of length. The length scale T may be used 
to define a second Peclet number P = vT/D which may 
be treated as a dimensionless measure of sample loading. 
A series of simulations are conducted with peak heights 
in the range 4> m = c m /c^ — 0.01 (low sample loading) 
to 0.8 (high sample loading). The initial co-ion concen- 
tration c p is assumed constant throughout the domain. 
Then the counter-ion concentration is determined by the 
local electro- neutrality constraint, Eq. ©. An infinite 
domain is approximated by a finite computational box 
of length much greater than wq. The values of the con- 
centrations are held fixed at the domain boundaries and 
d(j) e /dx is set to the constant value —E°°. The domain is 
chosen to be sufficiently large that the perturbations of 
the concentrations and fields are always negligible near 
the domain boundaries. 



A. Model System 

We will consider a three ion system consisting of sam- 
ple ions, co-ions and counter-ions. Results will be ex- 
pressed in terms of dimensionless variables: all lengths 
are in units of a characteristic length wq determined by 
the initial peak width, time is in units of wq/v, where v 
is the migration velocity of an isolated sample ion in the 
applied field (E°°) and the electric potential is in units of 
E°°Wf). All concentrations are in units of c£°, where is 
the concentration of negative ions in the background elec- 
trolyte. In order to define a minimal problem with the 
fewest possible parameters, we assume that the mobility 
(u) is the same for all the species, and therefore, so is the 
diffusivity (D), in accordance with the Einstein relation 



B. Numerical Method 

We will solve the governing equations for ion transport 
in solution, which are 



dcj 
dt 
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dx 



d(j) e 

dx 



Oct 
dx 



= 



(2) 



where Ci is the concentration of species i (i — 1, 2, . . . , N) 
with electrophoretic mobility /ii and diffusivity Di. 
Electro-osmotic flow is neglected so that the problem 
is one dimensional and may be described using the co- 
ordinate x along the capillary and time since injection, 
t. On account of the requirement of local electroneutral- 
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FIG. 1: Time evolution of the profiles of the normalized 
sample ion concentration ((/>), electric field (E/E<x>) and 
Kohlrauch function (K/Koo) in the case of weak sample load- 
ing (<f> <C (f) c )- The Kohlrausch function spreads only by dif- 
fusion so that the sample peak rapidly migrates into the zone 
where K — Koo (movie online). 
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(zi is the valence of the ith species). The electric po- 
tential (p e may be found from the equation of current 
conservation: 



0_ 

dx 



N a j. 

ZifAiCj, 

i=l 



dx 
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dcj 
dx 



0. (4) 



Eq. (j4| may be readily integrated to yield the local elec- 
tric field, E = —d x <j) e : 

£ (M)= £ °° £ ^f (5) 

where the superscript oo indicates the value of the respec- 
tive variable far away from the peak and the summation 
is over all species. 

A finite volume method is used to discretize equations 
@ and (HJ) in space using an adaptive grid refinement 
algorithm that is enabled by applying the Matlab li- 
brary "MatMOL" [l(J. The spatially discretized system 
of equations is then integrated in time using the Matlab 
solver "ode45" [TTI] which is based on an explicit Runge- 
Kutta (4,5) formula. Equations ^ and ((H) automatically 
ensure that the electro- neutrality condition, Eq. (O, is 
satisfied and this is verified at each time step. 



C. Results 

Figure 1(a) shows the profiles of the normalized sam- 
ple concentration <fi(x, t) = c n /d^ at fixed times vt/wQ = 



FIG. 2: The same as in Fig. 1 except here the amplitude 
of 4> exceeds cf> c - Here the part of the peak above the value 
<j> = (f) mteT appears to be effectively immobilized. The middle 
panel shows that the stagnant zone is due to a sharp reduction 
in the electric field caused by the very high electrical conduc- 
tivity in this zone. The assumption K = is clearly invalid 
as a part of the peak remains trapped in the injection zone 
(movie online). 



0, 0.5, 2.0, 4.0 and 8.0 for the case of low sample loading. 
Fig. 1(b) and (c) show respectively the profiles of the cor- 
responding electric field E(x, t) and the Kohlrausch regu- 
lating function K (x, t) = (c p + c n + c)/u. The Kohlrausch 
regulating function is a useful quantity for describing 
electrokinetic transport. If all ionic species have the same 
diffusivity, K(x,i) evolves as a passive scalar If ionic 
diffusivities are treated as zero, then K(x,t) is a con- 
served quantity [l2j |. It is seen that K(x,t) remains lo- 
calized near the injection zone and spreads only slowly by 
molecular diffusion. The sample peak on the other hand 
moves to the right and after a short time, the sample 
peak essentially lies in a zone where K = K^. This illus- 
trates the behavior postulated earlier that makes possible 
a simplified description in terms of the one dimensional 
nonlinear equation 
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d_ 

dx 



1 — acf> 



D 



dx 1 



(6) 



If <f) is small, Eq.© reduces to the Burgers' equation on 
Taylor expansion of (1 — acj))^ 1 . In the vicinity of the 
sample peak, the electric field is functionally related to 
the normalized sample concentration; E — E°° / (1 — ot(f>). 
Here a is the "velocity slope parameter" introduced in 

1- _ 

It may be shown [5j that the requirement of positivity 

of co- and counter-ion concentrations implies that only 
sample profiles satisfying the condition (f> < cf> c , where <p c 
is a positive number, may be described by the theory. 
We will call such profiles "realizable". The critical con- 
centration, 4> c , is given by <fi c — (z p — z n )/{z p — z) when 
z < and 4>c — —[z n (z p — z n )/[z p {z — z„)] when z > 0. 
When the parameter a > 0, it may be shown with some 
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FIG. 3: The normalized rate of change of variance as a func- 
tion of dimensionless time vt/wo for (a) Pe = 100 and (b) 
Pe = 200 for three different values of sample loading (P) . At 
long times, the peak is seen to spread with an effective diffu- 
sivity D eS = \im t ^o {2D)' 1 (da 2 /dt) given by Eq. (37) of [| 
and indicated here by the horizontal dashed line. 



FIG. 4: The normalized convergence time vt*/wo where t* 
is the time taken for the effective variance {2D)- 1 (da 2 /dt) 
to reach 95 percent of its final asymptotic value of D e g as a 
function of the sample loading (P) for different values of the 
Peclet number (Pe). 



simple algebra that <fi c < tfi' c = or 1 (see Appendix), so 
that the singularity implicit in Eq. (|6j) when <fi = <p' c is 
never reached for realizable solutions. Fig. 2 shows the 
behavior of the system for initial conditions that are not 
realizable. In this situation, a stationary "barrier" devel- 
ops at a fixed spatial location corresponding to a certain 
value (j> = (f> mter < (f) c . The sample ions move more or less 
freely on crossing the barrier but are effectively immobi- 
lized on the left of the barrier. This is due to the greatly 
reduced strength of the electric field in the injection zone 
where the electrical conductivity is high. This is clearly 
seen in Fig. 2(b) which shows a sharp reduction in the 
electric field in the injection zone. Only sample ions near 
the edges of the zone are able to "leak out" and are car- 
ried to the right as an advancing wave. Since part of the 
sample profile remains quasi stationary, the assumption 
of the constancy of the Kohlrausch function, K = 
can no longer be made for non-realizable concentrations. 
Thus, Eq.©, which would have led to unphysical neg- 
ative concentrations for such non-realizable profiles, is 
not applicable until after a sufficient time has evolved so 
that <f> is reduced to a value below (p mter throughout the 
domain. 

Fig. 3 shows the variation in time of the quantity 
(2D)~ 1 da 2 /dt for a series of different values of the dif- 
fusivity and sample loading characterized by the pair of 
Peclet numbers (Pe, P) . If the profile spread purely by 
molecular diffusivity, this quantity should approach one 
asymptotically. However, it is seen that the long time 
asymptotic value is not one but rather D & q which de- 
pends solely on P. The dashed line shows the theoretical 
value of -Dgff predicted by the weakly nonlinear theory 



based on solutions of the Burgers' equation [5]. Thus, 
once the system has evolved long enough, and dispersion 
has caused the amplitude to drop sufficiently, Burgers' 
equation provides a valid description of the peak evo- 
lution. However, a real separation happens in a finite 
capillary and the long time limit may not necessarily ap- 
ply. A quantity of interest is the timescale characterized 
by t„: the time needed for the quantity (2D)~ 1 da 2 / dt to 
relax to 0.95 of its asymptotic value -D e ff ■ If the separa- 
tion is conducted in a capillary of length L, the question 
of interest is whether t* is small or large compared to the 
total separation time T — L/v. In Fig. 4 we show the 
normalized time vt*/wo from a series of simulations with 
different values of (Pe, P). Clearly, vt*/wo is a monoton- 
ically decreasing function of P. This can be anticipated 
from the theory of nonlinear waves Q: the higher the 
amplitude, the quicker a shock or shock like structure is 
formed. In contrast to the effective diffusivity shown in 
Fig. 3 which depends on P but not on Pe, the time to 
reach the asymptotic state does depend on Pe. In fact, 
as Fig. 4 shows, the curve vt^/wo as a function of P is 
shifted upwards as Pe is increased. Indeed, larger Pe 
corresponds to lower diffusivity and therefore a longer 
time for the peak to spread and its amplitude to fall suf- 
ficiently for the weakly nonlinear description to be valid. 
Typical values of the physical parameters in a microchip 
based system may be wq ~ 100//m, L ~ 5 cm, so that 
vT/wo ~ 500. Thus, Fig. 4 suggests that the Burgers' 
solution does describe the peak dynamics for most of the 
separation time except for possibly a relatively short ini- 
tial transient. 
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FIG. 5: Schematic diagram describing approximately the ini- 
tial phase of peak evolution when (j> > <f> . The domain can be 
divided into an "Initial Zone" an "Interzone" and a "Back- 
ground Zone" . All the dependent variables are approximately 
constant within each zone but may undergo jumps across zone 
boundaries. 



D. Analysis 

An approximate theoretical determination of the con- 
centration (j) mter may be provided using the conserva- 
tion equations. The method of doing this is in fact en- 
tirely analogous to the "Moving Boundary Equations" 
(MBE) [13j for describing advancing fronts (e.g. in iso- 
tachophoresis), except, in this case, the front happens to 
be quasi stationary. The conceptual framework is illus- 
trated in Fig. 3. The domain is decomposed into three 
parts: the "Initial Zone" where the sample is injected, the 
"Background Zone" ahead of the advancing wave where 
all concentrations equal their initial values and an "Inter- 
zone" between them. All variables are assumed constant 
within each zone but undergo a discontinuous change 
across zone boundaries. The values of the variables in 
each zone are indicated in Fig. 3. The boundary between 
the Initial Zone and the Interzone is stationary whereas 
the boundary between the Interzone and the Background 
Zone moves to the right. The arrows indicate fluxes of 
ions across the stationary zone boundary. Conservation 
of these ionic fluxes require 

j^ini ^ini jointer ^inter 

flini^mi _ E inter gutter (g) 

where E represents the electric field and <f> represents 
the concentration (normalized by c£°). The superscript 
("ini" for the Initial Zone, "inter" for the Interzone and 
"oo" for the Background Zone) indicates the zone in 
which the variable is evaluated and the subscript {p for 
cation, n for anion and no subscript for the sample) iden- 
tifies the species. Therefore, 

prater = ^ini ^ini^nter (g) 

For the inter zone, 

jointer „oo / linter , linter , ±inter\ /., 

K = c„ [4> +(p p +(p n ) /u 



= K 00 =c^ ($? + 1) /u 

= c™(l-z n /z p )/u. (10) 

The electro-neutrality condition (valid in all zones) is: 

z p (j>p + z„4> n + z(j> = 0. (11) 

By combining Eq. (|9]) and (fTU|) and using the electro- 
neutrality condition we get an equation for determining 

imfer 

(1 - z/z p )<f> mter + r(l - z n /z p )cj> inter = 1 - z n /z p , (12) 

where the ratio (j)™ 1 / (f)" 11 — r is a, constant determined 
by the ionic composition of the injected zone. Solving 
the above linear equation for <fi mter we have 

In our numerical experiment the cation concentration 
was chosen to be uniform, so that 4> p m = </>^° = —z n /z p . 
Since z : z p : z n = 1 : 2 : — 1, r = 1/0'"" — z/z n = 1.5, 
and, cj) mter = 0.55. This value is indicated by the dashed 
line in Fig. 2(a). Clearly, it correctly describes the con- 
centration of sample in the Interzone. 

Thus, the theoretical description developed in [5] may 
be used in the Interzone (cf> < <f> mter ) but not in the 
Initial Zone. In order that all ion concentrations be non- 
negative in the Interzone we must have (f> mter < <\> c < <fi' c . 
This inequality is indeed true as can be shown by some 
simple algebra (see Appendix). 

III. CONCLUSIONS 

The development of nonlinear waves in capillary elec- 
trophoresis in the limit of low as well as high concen- 
tration of sample ions was studied by numerical integra- 
tion of the governing equations. An idealized minimal 
model was considered consisting of a three ion (sample, 
co- ion and counter-ion) system of strong electrolytes [2C| . 
This study complements an earlier paper by the authors. 
There it was shown that, in the weakly nonlinear limit, 
the evolution of the sample concentration may be reduced 
to Burgers' equation, which admits an exact analytical 
solution. 

Numerical simulation revealed that the evolution of 
the peak proceeds in a way that is qualitatively differ- 
ent when the sample concentration is high. As a con- 
sequence of the sharp reduction of the electric field in 
the region of sample injection, the ion migration velocity 
in this zone is very small. Ahead of this zone the ions 
form a surge front with a step-like profile propagating to 
the right. This state of affairs continues until the dimen- 
sionless ion concentration (0) in the injection zone drops 
sufficiently so that <j> < (f> mter . The subsequent dynamics 
then proceeds in accordance with the weakly nonlinear 
theory ||. The value of cj) mter may be approximately 
calculated by using a simple model based on conserva- 
tion of ionic fluxes. 
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This qualitative change in the dynamics of peak evolu- 
tion explains the breakdown of the weakly nonlinear the- 
ory when the concentration <f> exceeds the critical value 
4> c - When exceeds a certain value mter < <^ c part 
of the propagating wave is effectively immobilized in the 
injection zone. It is then no longer correct to assume |5[ 
that the sample pulse would quickly move out to a region 
where the Kohlrausch function is constant. 

The model studied here is clearly oversimplified. In 
particular, real electrophoresis buffers contain many 
more than three ions including one or more weak acids or 
bases to maintain a stable pH. Further, complex effects 
due to inhomogeneities in the electroosmotic flow may 
be relevant In this paper we ignore these complex- 
ities and attempt to produce a detailed understanding 
of a "minimal" model problem. One may question 
whether the strongly nonlinear regime considered here 
is of relevance to actual laboratory practice. The 
answer depends on the numerical values of the critical 
concentrations (f> mter < C < <fi' c . If the sample and 



carrier ions have similar valences then all of these critical 
concentrations are of order unity. Thus, to exceed these 
critical values the sample ions in the injected plug will 
need to be present at concentrations approaching that 
of the carrier electrolyte. Such high concentrations 
are rarely employed in laboratory practice. However, 
if the sample is a macro-ion the critical values may 
actually be quite small. For example, at pH 2.0 Bovine 
serum albumin has a valence, z ~ 55 [15|. Then, in a 
univalent carrier electrolyte we have <fi c ~ 0.04, so that 
the strongly nonlinear regime studied here may be easily 
reached. 
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Appendix: Proof of the inequality 



< 0c < 



The critical concentration 4> c is defined as || 



if z < 
if z > 



whereas 



1 

o 



(z - z n )(z - Zp)' 



(A.l) 



(A.2) 



We need to show that cf> c < <fi' c when a > 0, that is, when 
z p > z > z n . To do this, evaluate the ratio <j) c /<l>c when 
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z p -z n {z-z ){z-z ) _ -Zn+Z < 1 
z p — z z n (z p — z n ) — z n 

2/i Zp Z n {z Z n )(z -Sp) Zp Z 



if z < 
< 1 if z > 
(A.3) 

which completes the proof. 

To prove the remaining inequality, (j> mter < </> c , we first 
show that r > ~z/z n when z < 0. To do this, we use 
the electro-neutrality condition to express §™ % in terms 
of the other variables 



iini 



-(-z-rz n ). (A.4) 



Now we must have (j)™ 1 



> 0. This is always true if z < 0, 
but if z > then we require that r > —z/z n . 
First suppose that z < 0. Then 



< 



(z p - z)/(z p - z n ) 



Now suppose that z > 0. Then 



linter 



1 



< 



(A.5) 



r+(z p - z)/(z p - Zn) 
1 

~(z/z n ) + (Zp - Z)j (Zp - Z n ) 



Zp Z Z n 



(A.6) 



linter 



i 



r + (z P - z)/{z p - z„) 



Thus, in all cases, (p mter < C which completes the proof. 



